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ABSTRACT. Several algorithms are presented for the accurate computation of the leaves in 
the foliation of an ODE near a hyperbolic fixed point. They are variations of a contraction 
mapping method in [25] to compute inertial manifolds, which represents a particular leaf 
in the unstable foliation. Such a mapping is combined with one for the leaf in the stable 
foliation to compute the tracking initial condition for a given solution. The algorithms are 
demonstrated on the Kuramoto-Sivashinsky equation. 

Introduction 

The Hartman-Grobman Theorem provides a local foliation for an ODE near a hyper- 
bolic point; through each nearby (base) point there is a pair of leaves that define a conju- 
gacy to the linearized flow. In the classic case where the base point is the hyperbolic point 
itself, one leaf is its unstable manifold, the other its stable manifold. In that case the leaves 
are invariant; for a general base point they are not. They can, however, be characterized 
by the exponential growth/decay rates of the differences between solutions that start on 
them. If the gap in the spectrum of the linear part sufficiently dominates the Lipschitz 
constant for the nonlinear part in a large enough neighborhood, and the spectrum is posi- 
tioned properly, the unstable manifold is an inertial manifold. Each solution is attracted at 
an exponential rate to a particular "tracking" solution on the inertial manifold. We present 
several algorithms for the accurate computation of the leaves in the foliation and as well 
as for the tracking initial condition for a given solution. The algorithms are demonstrated 
on the Kuramoto-Sivashinsky equation, which is an amplitude model of thin film flow (see 
[13] and references therein). 

There has been considerable analysis of foliations in the literature. The finite dimen- 
sional case was studied in [17], followed by treatments for particular partial differential 
equations (PDEs) in [19] and [1]. The exponential tracking property of inertial manifolds 
was established in [11]. We consider here the general Banach space setting, as in [2], [3], 
and [5], and follow the particular framework in [4]. 

The computation of different elements in a foliation have been treated separately with 
a variety of approaches. The survey paper [18] discusses a great number of methods for 
classic stable and unstable 2D manifolds. Those manifolds are global, but are generally 
not the graphs of functions, unlike inertial manifolds, which are usually assumed to be 
both [10]. Approximate inertial manifolds (see e.g. [13]) for dissipative PDEs are explicit 
expressions for the enslavement of the high modes in terms of the low modes. To reduce 
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the error in their approximation one must increase the number of low modes, and hence 
the dimension of the manifold. Direct computation of a global inertial manifold of fixed 
dimension was carried out in [24], while accurate evaluation of the enslavement at individ- 
ual low mode inputs was achieved in [15, 21]. The first efforts to compute tracking initial 
conditions appear to be in [22, 7, 23], in the context of center manifolds. Those methods 
involve an expansion of normal forms and iterative procedure to project onto a local basis 
of the tangent space of the center manifold. 

We present here a unified approach to computing all these elements. It is based on the 
Lyapunov-Perron contraction mapping on spaces of functions in time used in [4], and out- 
lined in Section 1 . The fixed point of the mapping is a particular solution of the differential 
equation whose initial value provides an enslavement of either the high modes in terms 
of the low, in the case of a leaf in the unstable foliation, or vise-versa in the the case of 
the stable foliation. These leaves are manifolds and graphs of functions whose Lipschitz 
constants are less than one if the linear term in the equation sufficiently dominates the non- 
linear term. Both functions are then combined to form yet another contraction mapping 
whose fixed point is the intersection of the manifolds. In the particular case where the leaf 
in the unstable foliation passes through a steady state (as in the inertial manifold), the fixed 
point of the combined contraction mapping is the tracking initial condition for any point 
on the leaf in the stable foliation. 

The key then, is to discretize the Lyapunov-Perron contraction mapping. This was first 
done in [25] for the particular case of the inertial manifold by using piecewise constant 
functions over increasingly finer time intervals. We adapt that approach to the stable foli- 
ation which requires an inner integration of the differential equation, which happens to be 
forward in time, so it is practical for PDEs. This is done in Section 2. We then consider 
improvements based Aitken's acceleration and Simpson's method. The new methods are 
applied to the particular case of an inertial manifold in Section 3, and then combined for 
stable and unstable foliations to compute tracking initial conditions in Section 4. We wrap 
up in Section 5 with a comparison of long time dynamics of computed tracking initial con- 
ditions and a linearly projected initial conditions. Public domain software for this approach 
is available at [6]. 



We recall here the main features of foliation theory for ODEs, following the presentation 
by Castaneda and Rosa [4], where the proof can be found. Let X and Y be Banach spaces 
mdZ — X x Y be endowed with the norm ||z|| = || (x,y)\\ = max{\\x\\,\\y\\}. LetF : Z — >X 
and G : Z — > Y and A, B be two linear bounded operators defined on X and Y . Consider 
the following system of ODEs: 



1. Assumptions and Foliation Theory 



(1.1) 



i = Ax + F (x,y) 
y = By + G(x,y). 



The assumptions are the following: 



(1.2) 
(1-3) 
(1.4) 



j|e A '|| <e at and ||e" B '|| <e~ pt , Vf > 0, and for some a,j3 el. 
H(z) = {F{z),G{z)) is Lipschitz with Lip(H) < 8 and H{0) = 0. 
25 < j3 — a, called the spectral gap condition. 
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Since H(z) is Lipschitz and A, B are bounded operators, it is known that the autonomous 
differential equation 

(1.5) z = Cz + H{z), 

where C = A x B possesses a global unique solution for any given initial condition. We 
denote by z(t,zo), the solution of (1.5) with initial condition z(0) = zo G Z, andx(f ,zq) and 
y(t,zo), the X and Y components of z(t,zo), respectively. Thus, z(t,zo) = (x{t,zo),y{t,zo)), 
for all t G K and zo G Z. 

Typically, the nonlinear terms in most physical models are not globally Lipschitz. If 
this is the case and the system is dissipative, the nonlinear terms can be truncated outside 
the absorbing ball. More precisely, let p be the radius of the absorbing ball. Consider the 
prepared equation 

(1.6) z = Cz + H p {z), 

where H p :Z^Z which agrees with H for ||z|j < p and is globally Lipschitz. Since the all 
the long time behavior of the original system is in the absorbing ball, such a preparation 
leaves that behavior unchanged. One choice of H p is 

(1-7) H p = d p (\\z\\)H(z), 

with 

2 fl, fori G [0,1], 

(1.8) p (r) = 0(^ I ), 9(s) = I 2(i- 1) 3 -3(5-l) 2 +l, for s G [1,2], 
P [o, for i > 2. 

The main result in [4] is to characterize foliations by the exponential growth/decay of 
the difference of any two solutions with initial data in the same leaf. 

Theorem 1.1. (Foliation Theorem) 

(1 ) ( Stable) Z = \J ye y ~4(y, where 

(/) Jt y = {zo G Z :\\z(t,zo)-z(t, (0,y))\\ < \\z - (0,y)\\e^ +s \ V t > 0} 
(ii) . //., = graph(<P y ), for some <J> V : X — > Y such that ^(O) = y. 

(2) (Unstable) Z = \J xeX ^x, where 

(0 ^ x = {zoeZ:\\z(t,z )-z(t,(x,0))\\ < \\z - (x,0)|| e ^- 5 > V t < 0} 
(ii) jY x = graph^x), for some : Y -> X 5mc/i f/iaf ^(0) = x. 

(3) Both <t> v ant/ /tave Lipschitz constants bounded by 8/ (J5 — a — 8). 

The terminology stable foliation (unstable foliation) comes from the classic case where 
a < < /3, in which (jV§) are respectively the stable (unstable) manifolds of 0. The 
framework, however, also applies if a < /3 <0or0<a<p\ Regradless, JK§ and jV§ are 
both invariant; 

and similarly for ,Aq. As a consequence, given any initial data zo G Z, properties (i) in 
Theorem 1.1 define distinguished solutions in and 

Proposition 1.2. (Exponential Tracking) Given zq G Z, there exists a unique Zq G JVq 

\\z(t,zo)-z(t,zZ)\\ < e (a+5) '||zo-4ll 5 Vr>0. 

Definition 1.3. The solution z(-,Zn) is called the exponential tracking of z(-,zo) and zj 
is called the tracking initial condition of zo- 



4 



Y.-M. CHUNG AND M. S. JOLLY* 



By Theorem 1.1, the entire foliation is established. The intersection of a leaf from 
the stable foliation with one from the unstable foliation is a single element in Z; more 
precisely, for each x\ EX, y\ g Y, jtf yy n jV Xy is a single element of Z. This is proved in 
[4] by showing that 

Z:Z^Z,Z:(x,y)^(W Xl (y),0 yi (x)), 

has a unique fixed point, which is the intersection of two manifolds. In Section 4, we 
will implement an approximate of the £ map iteratively to compute the tracking initial 
condition. 
If 



(1-9) 



a + 8 < 0, and dim(F) < °° 



then jVq is an inertia! manifold, i.e. an exponentially attracting, finite dimensional, Lips- 
chitz manifold. Some of the key features of the foliation are illustrated in Figure 1(A). 



(A) Foliations for (1.13). 
^ := z(tQ,Zi) for (' = 0, and 1. 





Original Sequence — i— 




Aitoi sequent — "— 


\ 






\ 

\ 



j 


SIMP 


SIMPGS 


1 


0.2699E-1 


0.9353E-7 


2 


0.3128E-5 


0.4832E-6 


3 


0.2435E-5 


0.483 1E-6 


4 


0.4836E-6 


0.483 1E-6 


5 


0.4830E-6 


0.483 1E-6 


6 


0.483 1E-6 


0.483 1E-6 


7 


0.483 1E-6 


0.483 1E-6 



(B) Number of iterations versus (c) Number of iterations versus 
the absolute error. the absolute error. 



Figure 1 . (A) jV§ is the inertial manifold and Zq is the tracking initial 
condition for zo- ^ zo is the leaf in the stable foliation through zo and 
has property that for any two points zo, Zl on ./iC and any time t > 0, 
z(f,Zo) and z(t,Zi) lie on the same manifold ^ z ^ t ZQ y (B) Performances 
of PWCONST and PWCONST along with the Aitken's accelerationtest 
for (1.13) with parameters p = 10, Zo = (1> 1) an d x = 3, (C) Perfor- 
mances of SIMP and SIMPGS for the same settings as (B). 



1.1. An Example: Test Problem. In this section, we will give an example that will be 
used to demonstrate the algorithms we developed in this article. Consider the simplest 
system: 

d 



(1.10) 

dt \y J \y 

The foliation of the linear system (1.10) consists of vertical and horizontal lines. In order 
to obtain a nontrivial foliation, we will apply the transformation T = Ti o T\ , where 



(1.11) 
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By an elementary calculation we obtain 



(1.12) 



■y- 



PVi+r 



tan _1 (I+ ■ 



After the transformation T, the new ODE can be written as: 

/ - \tar- 

— x 



dx 
dt 



y—jtan (x) 



y — -tan (x) 



p(l + (y-j- } tan-i(x))2y/2 p(l + (y- Ltan-^x)) 2 )^ 



(1.13) 



dy 
dt 



=y tan 

P 



-j^tan 



Han (x) 



P (i + {y-\tan-^x)yyn 



'(*) 



p( l + {y-fjta 



-l 



W)2)V2 J 



Note that p is a parameter that controls the Lipschitz constant. In fact, we can characterize 
the complete foliation. However, for the purpose of this article, we are interested in only 
the inertial manifold (invariant unstable manifold) and a leaf in the stable foliation. 

First consider the inertial manifold for the new system. Since the invariant unstable 
manifold for the original system is the y-axis, its image under the map T is the invariant 
unstable manifold for the new system. That is 



(1.14) 



= T 



v + - tan 

J p 



'l+v 2 



-1 



Since the leaf in the stable foliation through zo for the original system is the horizontal 
line, y = % for a given yo, the transversal manifold for the new system is 

l 



= T 



(1.15) 

By substitution, one obtains 



- tan 



y = yo + -tan ! (x). 



Given any initial condition fo = {x~o,y~o), the tracking initial condition for the original 
system is (0,/o)- Thus, the tracking initial condition for the new system is 



(1.16) 



yo - 



±tan-!(- 



2. Computation of the Stable Foliation 

2.1. PWCONST algorithm. In [25], an algorithm is developed for the accurate compu- 
tation of inertial manifolds under the additional assumption (1.9). The main idea is to find 
the fixed point of the contraction mapping 

(2.1) W(y,y)(t) =e tB y+ f e^ A F(y(s))ds - f° ' e^ B G{\tf{sj) ds ,V t <0 , 
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on the Banach space 

(2.2) Sf ff = {vr6C((-oo,0],Z); ||y||<7 = sup e CT '|| v/(f)|| < «}, a 6 (a + 5,/3 -5). 

r<0 

The inertial manifold is the graph of the function: 

(2.3) W(y)=Pyr(0), 

where \j/ is the fixed point of and P is a projector from Z onto X. In the /* iteration, 
y is approximated by yW, a function that is piecewise constant on Nj time intervals of 
length hj. It is shown in [25] that \jf^ — s- \j/ as j : — >• °° provided iV/Ay — ^ 00 (e.g. = 2~ 7 
and Nj = j2->). In this approach the integrals in °i/ can be evaluated explicitly. 

To compute the leaf in the stable foliation through zo, we follow the existence proof in 
[4] and approximate the fixed point of the mapping 

(2.4) &o{<p,x)(t) = e tA x + J^e^ A [F((p(s)+z(s,zo))-F(z(s,zo))} ds 

(2.5) - j~ e^ B [G{<p{s) +z(s,zo)) - G(z(s,zo))} ds, 
on 

(2.6) ^ = {(peC([0,H,Z); !|(p|| (J = sup e - CT '||(p(r)|| <«.}, ae(a + 8,p-8). 

Let q> be the fixed point of The leaf in the stable foliation through zo is the graph of the 
function: 

(2.7) =yo + fi9(*-*o)(o), 

where Q = I — P. We modify the algorithm in [25] to fit & . The two main differences 
are: 1 . one needs to solve for the ODE forward in time — we use a 4-th order Runge-Kutta 
method (RK4); 2. an additional function evaluation is needed. The rate of convergence is 
linear so that we can use the Aitken acceleration process, as discussed in the next section, to 
gain a better approximation. Since piecewise constant functions are used in this algorithm, 
we denote it by PWCONST. 

2.1.1. Aitken Acceleration. Aitken's acceleration, also known as Aitken's A 2 process, is 
used for accelerating the rate of convergence of a sequence. The method works if one has 
a linear rate of convergence sequence. 

Definition 2.1. Given a sequence {-K/t}" =0 in K, Aitken's acceleration sequence, x^ ={) , 
is defined as 

(Ax*) 2 

(2.8) ^x k =x k -^-, 

where Ax^ := x k+ i —x k and A 2 X£ := Ax k + \ — &Xk- 

This is the classic Aitken acceleration for a sequence in R. We will also apply a vector 
version found in [20]. 

Definition 2.2. Given a sequence {Zjt}^ =0 in R", define 

Aza- := (z* + i -Zk,-,Zk+ n -Zk+n-i) and A 2 ^ := Az k+ i — Az k . 
Aitken's acceleration sequence for R" is defined as follows: 

(2.9) sfzk =z k - (Az k )(A 2 z k )-\z k+1 -z k ). 
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Note that both Azk and A 2 Zk are matrices of size n x n. 

The advantage of the Aitken sequence is that it converges much faster to the limit than 
the original sequence does. Moreover, computing the Aitken sequence is much cheaper 
than computing the original sequence because it is applied to elements in phase space 
Zk £ Z, rather than (a discretized version of) the function space on which ' , S? act. It 
amounts to post-processing the original sequence. A disadvantage is that it requires n+k 
iterations of the original mapping °l/ ( 3") to produce k terms in the Aitken sequence, which 
may be prohibitive if n, the dimension of Z, is large. We will see these numerical results in 
the next section. 

2.2. SIMP Algorithm. First, we give a recursive relation for 3 ZQ ((p,x) in the time variable 
. Second, we introduce the main algorithm. At the end of this section, numerical results 
for the test problem 1.13 will be given. 

Let {ti}^ be the uniform partition of a time interval and let h := f, + i — r,-. Let P be the 
projector from Z onto X, i.e. Pz = x, where z = (x,y) and Q — I — P. 

2.2.1. A Recursive Relation. We state the recursive relations in the following proposition. 

Proposition 2.3. 

( i) For z = 1, 2, . . . , N — I, the X-component of (<p, x) is 

P3 Z0 (( P ,x)(t I+l )=e^- l ^ A P.X ((p,x)(t i ^) + 

'' +l e^-^ A [F((p(s)+z(s,zo))~F(z(s,z ))}ds, 

where P ^^((p ,x)(to) — xand 

P.% (<p,x){ h )=e^ A x+ /" e^ A [F((p(s)+z(s,ZQ))-F(z(s 7 zo))} ds. 
Jo 

( ii) For i = N — 1, N — 2, . . . , 1, the Y -component of ,% a (<p,x) is 

e^ ( ( p,x)(f,_ 1 )= e ('-i-''+') B e^ ((p,x)(f, +1 )- 

ek-i-^[G((p(*) +z(s,z )) - G(z(s,zo))} ds, 

J ti-\ 

where Q^(<p,x){t N ) = - $~ e^-^ B [G{(p(s) +z{s,zq)) - G(z(s,zo))} ds and 
Q^(p,x){tN-i) = Q^,{(p^){t N )e^-^ B - f N e^-^ B [G((p(s)+z(s,zo))-G(z(s,zo))}ds. 



The proof of Proposition 2.3 is similar to the derivation of the recursive algorithm in 
[14]. 

2.2.2. The SIMP Algorithms. The integrals appearing in the recursive relation can be bet- 
ter approximated by Simpson's rule. In order to have a consistent order, one needs to 
adjust the first two approximations, i.e. P£? ZQ ((p,x)(to), P^K ((p,x)(t\), Q3K ((p,x)(tN), 
Q^{<p,x){t N -i). 

Since P^ ((p,x)(fo) = x, no error is introduced at t = tQ. When t =t\, 

P3 (<p,x)( tl ) = e'^ A x + J\^ A [F((p(s) +zM) ~ F(z(s,zo))}ds, 
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which involves the integral we need to approximate so that the error is the same as the one 
introduced by Simpson's rule. One can rewrite the integral as 

f 1 e^ A [F((p(s) +z(s,zo))-F(z(s,zq))] ds 

(2.10) = [ t3 e^ A [F((p{ S )+z(s,zo))-F(z(s,zo))]d S 

(2.11) - f ,3 e ill - s)A [F((p( S )+z(s,ZQ))~F(z(s,zo))]ds 

Jt\ 

For the integral (2.10), Simpson's 3/8 rule is used. It has the same order error as the classic 
form of Simpson's rule, which is used for the integral (2.1 1). 

The term Q,% (<p,x) (fjy) is the tail of the convergent improper integral. We will give an 
estimate on tjy so that the truncation error will be the same as the error of Simpson's rule. 
To do so, we need an estimate on <p, the fixed point of 3K ((p,x). 

Lemma 2.4. Let <p be the fixed point of £? zo . Then \\ (p\\ a < jzr^, where K = max{-^^, ^~^}- 
Proof. From [4], we have 

(2.12) ||^ (<p)||a<NI + K!l<Plla, 

for any <p £ Since <p is the fixed point of ,% , the left hand side of (2.12) can be 
replaced by ||<p||cj. By the gap condition (1.4), K< 1. Thus, by direct calculation, one has 



n / l|A ' 



l-K 



□ 



We are ready to give an estimate for <2^ (<p,A:)(fjv)- Note that 

\m zo (cp,x)(t N )\\ = || re^ B [G(cp(s)+z(s,zo))-G(z(s,zo))}ds\\ 
Jt N 

< I™ e {tN - s) P8\<p{s)\ds< e {lN - s) He as \\q>\\ a ds 
Jt N Jt N 

= 6\\<p\\ a M fe^-^ds = -^—e^Me 

p - a l-K" 1 — k 
Let h = t\ — to. Since the error for Simpson's rule is of order 0(h 5 ), we take 

-^--°*||jc|| <h 5 . 



l-K 

Hence, we choose fjy so that 

(2.13) e a,N < ^--^h 5 . 

k\\*\\ 

For Q£? Z0 ((p,x)(tN-\), the technique is similar to the one used for P Z?L{(p ,x){t\) . Rewrite 
Q37- ((p,x)(tN-i) as two integrals. 

«ftv-i-»)*[G(9(j) +z(s,zo)) - G(z( S ,z ))]ds 

r JV-l 
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(2.14) = e^-'-^ B [G((p( S )+z(s 7 zo))-G(z(s 7 zo))}ds 

Jt N _ 3 

(2.15) - f N ~ l e^-^ B [G((p(s)+z(s,zo))-G(z(s,zo))]ds 

J'N-3 

As before, the integral (2.14) can be approximated by Simpson's 3/8 rule; the integral 
(2.15) can be approximated by Simpson rule. Hence, the error for Q^ Z0 ((p,x)(tN-i) is of 
order <&(h 5 ). 

Before we present the algorithm, we remark on successive iteration. We have to ensure 
that the initial guess, q>°, is in the space & a defined in 2.6. In the PWCONST algorithm, 

is simply a constant function, which is in the space $# CT (defined in 2.2) if we may take 
0~ < (i.e. if a + 8 < < j3 — 8). However, in the case of the stable foliation, if (p° is 
constant, it may not be in J? a . Since ||<p°||c7 = sup f >oe _OT ||<p (f)||, ||(p°|| CT is not bounded 
if <7 is negative. Instead of constant <p°, we take 

cp°(x,t)=e at (x-xo). 

Then 

||fl> ||a = supe- a < || <p°(f) II = supe^-^'l^-jcoll < - 
r>0 t>0 

since (a — a) is negative by the gap condition (1.4), where a <G (a + 8, j3 — 8). 



SIMP algorithm 

Input: zo = (xo,yo), x, h, the step size, and J, the number of iterations. 
Output: 4> zo (x) 
1: Choose by (2.13) and t\ = ; x h for i = 0, , 1 , N. 

2: p° (*,/;) = (x-Xo)e at >. 

3: Compute the solution of the ODE, z(f;,zo) f° r i = 0, , 1 , Nby some ode solver. 
4: Evaluate Fj := F(z(ti,zo)) and G; := G(z(ti,zo)) for i = 0, ,1, N. 
5: for j = -> / - 1 do 

6: ^:= J F(^'(f,)+z(f„zo))-^,for/ = 0, ,1, ,...,iV 
7: G/ := G((pj(ti)+z(ti,zo)) - G u for i = 0, , 1, , N 

8: P& Zo ((pj+ 1 ,x)(t )=X-X0 
9: P^ Z0 (c P J+\x)(h) = 

e tlA (x~x Q )+S38(e^- , ^ A F l l;k^0,l,2,3)~S(e^- , ^ A F^,k^ 1,2,3) 

10: e^zo(<P 7+1 ^)(^)=0 
11: Q% W +1 ,x)(t N -i) = 

-S38(e^-^ A G' k ;k = N -3,N -2,N - l,N)+ 

S(e ( fN-^- t k) A G{;k = N-3,N-2,N- 1) 
12: for i = 1 -> A' 1 do 

13: P^ (^ +1 ,x)( f! - +1 ) = 

e( ( >'+i-'>--i) A J D ^ (9^',x)(f I -_i) +S(e ( ?>+ l -ri A Fj;k = 

14: e^ (^+ 1 ,x)(f w _,_ 1 ) = 

e ^- i -i-<A r -m)A Q ^( ( p; )X )( % _. +1 )_ 

5( e (%-,-i-^M G >. k = N — i—l,N—i,N — i + l) 
15: end for 
16: end for 

17: 4> ;o (x)=yo + e^ ((p J ,x)(0) 
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For the notations in the SIMP algorithm, 

S(f k ; k=i,i + l,i + 2):=^ (fi + 4f m + f i+2 ), 
S38(f k ;k = i,i+l,i + 2,i + 3):= j^(/< + 3/ I+1 + 3f i+2 +/,+ 3 ), 
where h = f, + i — ?,-. In particular for Simpson's rule in the /-loop, we have 

5 ( e fe+l -'k)A F j. k = j _ 1, j, j + 1) 

= ^ e fe + ,-^l)A (F( ^ (fi _ l)+z(fi _ hZo)) _ F(z(f ._ 1;Zo))) + 
4e (', + 1 ^)A (F((p 7 (r . ) +z ( thZ0 )) -F(z(t U Z0)))+ 

e (^- f .-+i)A (F(<p y (t . +l)+ ^. +1)Zo)) _ j p (z(? . +1)Zo )) )] . 

Since f^ (^) 7+1 ,jt:)(f,_i) and P^ (<P"' +1 , jc) (f,) are computed in the previous steps, we 
could, as in Gauss-Seidel iteration, use these to obtain a better approximation. More pre- 
cisely, we do the following: 

-<k)A F L k = i _ lii+l ^ 

Ae {t ^ A {F(P^ + \ ti ) +x(t i ,zo),Q<P J (t l ) +y(ti,zo))-F(z(t u zo)))+ 
e (n + i-ti + i)A^ F{9 j {tj+l)+z(ti+uZo)) -F(z(t t+1 ,zo)))]. 

A similar partial update can be done in reverse for 

S (e fa-'- 1 -' k )A G{ ; k = N - i - 1 , N - i , N - i + 1 ) . 
We call the resulting algorithm SIMPGS. 

2.3. Numerical results for the Test problem. We apply PWCONST, PWCONST+Aitken, 
and SIMP to compute a leaf in the stable foliation for the test problem with p = 10. We 
take as the inputs for the algorithm 

(^o)=ni,i)=d +T ^ !J i+^- 1 (i +I ^ ! )) 



and x = 3 + -y^yj- The output should be 



>'=l + : 
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In Figure 1(B), the error of the original sequence decreases roughly by a factor of 1 /2 
as j increases; the error of the Aitken's sequence decreases roughly by a factor of 10 . In 
Figure 1(C), observe that first, both algorithms converge and both errors are of the same 
order. Comparing the two columns, there is a huge difference in the first two iterations. In 
1 or 2 iterations, the errors from SIMPGS seem to be saturated while errors from SIMP are 
saturated in 5 or 6 iterations. This is typical. Different h and different inputs give similar 
results. 
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3. Improved Computation of Inertial Manifolds 

Both algorithms discussed in the previous section can be adapted to compute an inertial 
manifold. In this section, we will show by the numerical evidence that these are the im- 
proved methods and we will apply those methods to the Kuramoto-Sivashinsky equation 
(KSE). 

The KSE with periodic and odd boundary conditions can be written 
du <9 4 m r d 2 u <9m, 

0.1) Tt +A w +Y[ W +u M ]=Q ' 

(3.2) u(t,$) = u(t,$+2n), u(t,-$)=-u(t,$). 

The solutions may be represented by the Fourier sine series 

u(t,$)= £ Uj (t)e^ = Y,bj{t)sin(j$), 

/=— j=i 

where the reality and oddness conditions on u give 

ib j 

ii-j = u j and Uj = — — . 

For the tests in this section, we use a 16 mode Galerkin approximation and take dim(X) = 
djra(T) = 8andy=32. 

It has been shown that KSE (3.1) has an inertial manifold (see [9]) and its lowest di- 
mension has been studied in [9], [26], and [16]. In particular at y = 32, it is shown in [16] 
that the computed global attractor is contained in a ball of radius 15 (in the L 2 -norm) which 
when used in the preparation in (1.7) yields an inertial manifold of dimension five. In what 
follows we take dim{X) = dim(Y) = 8, as the larger gap at this splittinbg leads to more 
rapid convergence for all methods. Since a limit cycle is contained in the global attractor, 
which in turn is on the inertial manifold, we pick a test point, mo = xo +yo, on a limit cycle 
and pick the low mode component, yo, as an input of the algorithm and test how well we 
recover the high modes, xq. 

3.1. PWCONST plus Aitken's Acceleration. As we mentioned before, since the rate of 
convergence is linear in PWCONST, we can apply Aitken's A 2 process to accelerate the 
convergence. The reason there are only 5 iterations in the Aitken's sequence in Figure 
2(A) is that in calculating one iteration in the Aitken's sequence, one needs n + 1 iterations 
from the original sequence, where n is the dimension of the elements of the sequence. The 
error of Aitken's sequence is much better than the original sequence. Another advantage 
of Aitken's A 2 process is the cheap computation. In order to obtain the similar error with 
the original sequence, one needs to compute more iterations, which is more expensive than 
computing the Aitken's sequence. To observe the difference between these computational 
efforts, we compare the number of multiplications for the two algorithms, namely PW- 
CONST and the Aitken's sequence ignoring multiplications needed for the evolution of 
nonlinear terms and the integrals. For the PWCONST algorithm, the required number of 
multiplication is at least 

l/ j=1 j x 2' x 6 x dim(Z) . 

In this case, dim(Z) = 16. If we compute 10 terms in the PWCONST algorithm, the 
minimum number of multiplications is 

L)ij x2 j x6x dim(Z) = 1769664 . 
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i 

j 


xk 







i 


0.616E-4 


N/A 


2 


0.321E-4 


N/A 


3 


0.146E-4 


N/A 


4 


0.707E-5 


N/A 


5 


0.348E-5 


N/A 


6 


0.173E-5 


N/A 


7 


0.862E-6 


N/A 


8 


0.430E-6 


N/A 


9 


0.214E-6 


1.595E-11 


10 


0.107E-6 


1.629E-11 


11 


0.537E-7 


1.627E-11 


12 


0.268E-7 


1.627E-11 


13 


0.134E-7 


1.627E-11 



(A) Performances of PWCONST 
and PWCONST + Aitken's accel- 
eration. 




2 



(B) Performance of the SIMP algorithm. 



FIGURE 2. Improved computation of the inertial manifold for the 16 
mode Galerkin approximation of KSE with y = 32. (A) Absolute error 
table for two algorithms. x J (stx^) is generated by PWCONST (PW- 
CONST+Aitken). (B) We fix (h,N) = (le - 6, 100000) and vary dim{Y). 



On the other hand, to compute one term of the Aitken's sequence, one needs to compute 
dim(X) terms in the PWCONST algorithm, solve a linear system, and calculate a matrix- 
vector multiplication. Thus, the number of multiplications for computing two terms of the 
Aitken's sequence is 

dim(X) + \ 

(3.3) £ jx2 J 'x6x dim(Z) + (dim(X) 3 + dim(X) 2 ) x 2 = 787776. 

j=i 

3.2. The SIMP Algorithm. To compute the inertial manifold and a leaf through zq in the 
stable foliation is to compute the fixed point of in (2.1) and ,% in (2.4) on the Banach 
space & a in (2.2) and in (2.6) respectively. These two maps are similar and in fact, 
if we formally replace (zo,(p,t,x,A,B,F,G) in the 3~ map by (Q,\j/,—t,y,B,A,G,F) and 
drop the z(s,Zo) terms, we obtain the "i/ map. With this observation, we can easily modify 
the previous algorithms to compute the fixed point of the % map. 

There are two parameters in SIMP, namely step size, h, and the number of points, AT. 
Given h, one can choose Af by (2.13), though this requires an estimate on the Lipschitz 
constant of the nonlinear term, such can be found in [16]. Here, we will choose h and 
experimentally. 

Next, we investigate how the spectral gap affects the algorithm. For the KSE, since 
the linear term is a diagonal matrix, the gap is the difference between two consecutive 
eigenvalues. 

When dim(Y) = 1 and dim{Y) = 2, the sequences are convergent but they are not con- 
vergent to the inertial manifold. When dim(Y) = 3, the convergence is slow; it took about 
50 iterations to reach the saturated error, which is about 10~ 8 . The reason could be that the 
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gap condition is barely satisfied for the region of phase space visited by the algorithm. As 
we increase dim(F) and hence the gap, the convergence is faster as Figure 2(B) shows. 

4. Computation of Tracking Initial Conditions 

4. 1 . Algorithm for tracking initial condition. It is shown in [4] that the exact tracking 
initial condition of a base point (jci ,y\ ) is the fixed point for the mapping 

(4.1) E: (x,y)^(V Xl (y),<t> yi (x)). 
We will fix j\ and j 2 and iterate the map 

(4.2) J\7W2 . (^y) _> i^h^ (y^h 

As is the case for E, E ;i j2 is a contraction mapping and the fixed point, (x* ,y*), is the 
intersection of the two manifolds that are the graphs of *¥ }l X \ ar, d *^yi ( x )- 

4.2. Convergence of the algorithm. The following results show that the algorithm for 
the tracking initial condition converges under a stronger gap condition. 

Lemma 4.1. ||(p y (xi) - (p j (x 2 )\\cr < t^II x i ~ x 2||, where K = max{j^, ^r^} and 
||<P||o- := su Pj .> e— CTf ||<p(r)||. 

Proof. For j > 1, define <p 7 recursively by 

<p-'(x)(f) = 3r^-\x)(t), where <p°(x)(t) = e at (x-xo), Vf e [0,«). 

Now, let f S [0,o°) and consider 

H^CxOCO-^'C^coil = ||^ (^'- 1 ^i)(0-^o(<p^ 1 ^2)(?)ll 

=\\e tA (x 1 -x 2 ) + r e ('-^[ J F(^- 1 (x 1 )W+z(^Zo))-^(^'- 1 (^)W+^,Zo))]^ 
Jo 

> e ('-^[G(^- 1 (x 1 )(5)+z(5,zo))-G((p^ 1 (x 2 )( i )+z( i ,zo))]^|| 

<m fl jc{ e ffi | |*i - x 2 \ | + f e^ a 5 1 1 q>j~ 1 (x i ) (s) - <p j - 1 (x 2 ) (5) 1 1 ds, 
Jo 



Let (7 € (oc + — 5) and multiply the last inequality above by e ct to obtain: 

e-^llp'faXO-p'teXOII 

< maxie^-^'Wxi -x 2 \\ +e-°< f e^- s)cc 8\\(p J - l { Xl )(s) - (p J - 1 (x 2 ){s)\\ds, 

Jo 

e-°' jf° e ('-'>/» 5 1 1 1 (jc 1 ) (s ) - - 1 (x 2 ) ( j) 1 1 d*} 
^mflr^-^lljti-jell + SII^-'fxiJ-^feJIla [' e {a - a ^'- s) ds, 

JO 

^lljci-xall + icll^-H*!)-^" 1 ^)!^ 
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Now, since the right hand side of the inequality does not depend on f, take the supremum 
over t and obtain 

- 9 J(x 2 )\\a < -x 2 || + K\\<P J -\xi) - <P ] -\x 2 )\\ a 

< ||*i -X 2 \\ + K\\X! -X 2 \\ + K 2 \\(p j - 2 (xi) - (p i ' 2 {x 2 )\\ (J 

< \\xi -X 2 \\ + K\\ Xl -X 2 \\ +... + K y ||p°(jCl) - (p Q {x 2 )\\a 

= (l + K-+K- 2 + ... + ^')!|x 1 -X 2 || < ; 

□ 

Since 4>- (x) = yo + Q(p(Q), we can obtain an estimate for 4> ; . For the exact manifold, 
from [4] we have Lip(<i> zo ) < ^ _ 5 a _ s - 

Lemma 4.2. 

W{xi)-&{x 2 )\\< p 8 _ al l _ K \\ Xl -x 2 \\. 
Proof. Apply Lemma 4. 1 to find that 

ll^xi)-*^)!! = || re- sB [G(( P j-\x l )(s)+z(s,z ))-G(( P j-\x 2 )(s)+z(s,z ))]ds\\ 
Jo 

< r e-'t S\\ (pj- 1 (*i) (s) - (p^ (* 2 )(*)| I ds 
Jo 

= f e-^e sa 8e- sa W-\x x ){s) - cp^\x 2 ){s)\\ds 
Jo 

<8W-\ Xi )-^-\x 2 )\\ c fe^-^ds 

Jo 

= ^W-\^)-^-\x 2 )\\ a < ^ 8 _ al \ K \\ Xl -x 2 \\. 

□ 

We have an analogous result for *F, the proof of which is similar. 
Lemma 4.3. 

||yy(y l) _yJ C y 2) ||< _|__l. || yi _y 2 ||. 

Since we now have both estimates for <t> ; and ^P 7 , ready to show that J2 is a contrac- 
tion mapping. 

Proposition 4.4. If 

(4.3) 48<l3-a, 
Y. JlJ2 defined above is a contraction mapping with 

Lip(lZ h ' h ) < — —, 

1 — K 

Proof. By Lemma 4.2 and 4.3, 

Hs'^Cjci.yO-^'^C^ya)!! = moxiii^OcO-^C^ii.ii'F-'CyO-^Cya)!!} 

< -, max{\\xi -x 2 \\,\\yi ~y 2 \\}- 
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The condition < 1 is equivalent to K = max{-p^, ^z^} < 1 /2, i.e., 25 < O" — a and 
28 < p — a. The last two inequalities are equivalent to (4.3) for (7 = ^j^. □ □ 

Proposition 4.5. E- 71 '-' 2 has as its fixed point the intersection of the graphs of and<i> 12 . 
The proof is as for E in Castaneda Rosa [4]. 

4.3. Error estimate. Let Zq be the fixed point of EJ 1,j2 , and Zq the exact tracking initial 
condition (fixed point of E). We seek an estimate for ||zq — We start with two lemmas. 

Lemma 4.6. Let zq G Z. For any x G X and positive integer j, we have 

\\( P (x)-(p\x)\\ (J <^\\ (p \x)- (p (x)\\ (J , 
where (p is the fixed point of 3? zo ((p,x) and K = max{-^-^, ^~a}- 

Proof. Let m be an integer such that m > j. From [4], we know that Lip (3?) < K and by 
the gap condition, K < 1 . 

\\cp m (x) - 9>(x)|| a < ^=; + lll<P*M - ^ Wlla 

<2? =J+ llC t - 1 ||9 1 W-9 Wlla< T ^||9 1 W-9 Wlla. 
Since the last inequality is independent of m, take m — > °° to obtain 

IIVW-^WIIa^Y^II^W-^WIIa. 

□ 

Since 4>- (x) := yo + 6p(*)(0), the estimate for H^x) — < t >; (x)|| is straightforward. 
Lemma 4.7. Under the same assumption above, one has 

||4>(x)-^(x)||< T ^||( P 1 (x)-(p ( V)lla. 
Proof. By the definition of <t> and the previous lemma, 

||*(x)-*>(x)|| < ||fi9W(0)-j2^W(0)|| 

<\\ ( p(x)-^(x)\\ a <^\\ (p 1 (x)- (p °(x)\\ a . 

□ 

We have an analogous result for *F, the proof of which is similar. 
Lemma 4.8. Under the same assumption above, one has 

\\^{y)-^ ] iy)\\<^ c \W\y)-w\y)\\o. 

By Lemma 4.7 and 4.8, one can deduce that 
(4.4) 

||E(z)-E^(z)|| < ^max^Mlp'M-^MIk K^\\ W l {y) - W°(y)\\ a }, 
where j = min{j*i,/2}. We are ready to give an estimate for \\zq — Zq||. 
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Proposition 4.9. Let z$ be the fixed point o/E and be the fixed point ofV. Then 

K j 

114-411 - c >i^Y^2^' 
where j = mm{j l , j 2 } and c h , h = max{ ~i | | (p 1 (x) - p° (x) \ \ a , ~ j \ | y/ 1 (y ) - y/° (y) || a }. 
Proof. By (4.4) and Proposition 4.4, 

11^ — ^Sll = ll^(^) — ^'-^(^B)!! < ll^(^) — ^"-^C^)!! -k ||^'-^(^) — 



Therefore, we have 



\4-4\\<c hl j2 



1 — 2k 

□ 

For simplicity, we denote E' := E^ 7 in the following numerical tests. 
4.4. Application. 

4.4. 1 . Test Problem. In this section, we demonstrate the computation of the tracking initial 
condition for the test problem (1.13). In this test, (xo,yo) = ( 1 , 1) and p = 1 and hence by 
(1.16), the exact tracking initial condition is 



(4.5) 



«oW life \„ (7.0711* -02 



10 La " V 10v /2 

By Proposition 4.4, for fixed j, the Lipschitz constant for E 7 ' is 



K 



l-K 

In the test problem, a = — 1, J3 = 1, and 8 < l/p and hence fc < I /p. In this case, = 10, 
and therefore Lip{L J ) < 1 /9. As Table 3 (A) shows, in the first two iterations we observe 
that the error decreases roughly by a factor of 0.01 which confirms Lip(U) < 1 /9. 

To verify Proposition 4.9, we will vary j and observe how the error depends on j. As 
in Table 3 (B), we can see that the error decreases roughly by a factor of 0.1 as j increases 
which also confirms Proposition 4.9 that the rate of the convergence is approximately 
K(«l/10). 

4.4.2. Approximate Inertial Form of the KSE. We express the KSE in the functional form: 

— +Lu+R(u) = 0, u G Jrff , 
dt 

where Jff is an appropriate Hilbert space (see [13]). The linear operatorL is given by 

along with periodic, odd boundary conditions. The remaining terms are then collected in 
R. The infinite-dimensional phase space is split into low- and high-wavenumber modes 
by means of the projectors 

P =P n : Jt? —> span{sin(%) ,sin(2%) , . . .,sin(nt;)}, Q = Q n =I — P n - 
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3.5429 


E 1 


4.97E-4 


z 1 


2.4E-2 


E 2 


8.49E-6 


z 1 


8.4E-5 


E J 


1.85E-7 


z J 


2.9E-7 


E 4 


4.50E-9 


z 4 


1.0E-9 


E 3 


1.16E-10 


z s 


3.3E-11 


E b 


3.25E-11 
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(A) 



(B) 



(C) 



Figure 3 . (A) This table shows the convergence of the tracking initial 
condition for (1.13) using SIMP for fixed j = 15 and z i+1 := E^(z' ) with 
h = 0.01. (B) This table shows the relation between j and the absolute 
error. (C) Convergence for data in (A) in phase space. 



Thus, u = p + q, where p = Pu and q = Qu. We will use the approximate inertial manifold 
(see [8]), <t>i (p) = —Lr x QR(p) . The approximate inertial form is 

(4.6) ^E-+Lp + PR(p + <t> l (p))=0, 

at 

which in our foliation framework would meanz = p,C = —L, and H(z) = —PR(z + <i>i (z)). 
Compare to the Galerkin approximation which amounts to replacing <t>i with <t>o = 0. We 
fix n = 3 as it was demonstrated in [13] that this is sufficient for (4.6) to capture the long 
time dynamics of the KSE for y g [0,36]. There are two reasons convergence should be 
slower for this reduced system. First, the gap in the eigenvalues is smaller. Second, the 
composition of the nonlinear term with itself in (4.6) makes for a larger Lipschitz constant. 

Given two initial conditions for (4.6), denoted by z\ and z%, which are close to each other 
but on opposite sides of a separatrix, we will construct the leaves of the stable foliation 
through each point. Since each leaf is an equivalence class, the points on the same leaf 
should have the same long time behavior. This property is demonstrated in Figure 4. 

5. Comparison between the tracking initial condition and the projected 

one 

Given an initial condition zo, it is natural to take zq := (0,yo), the linear projection onto 
the unstable eigenspace, as the approximate tracking initial condition since Y is tangent to 
the inertial manifold. However, the projected initial condition may have the wrong long 
time behavior as Figure 5(A) shows. Moreover, even if it captures the correct long time 
behaviors, the tracking rate is not optimal compared to the tracking rate for the tracking 
initial condition. By Proposition 1.2, one has that 

(5.1) ||z(f,zo)-z(f,z + )||< e (a+5), ||x -x+||,Vf>0. 

Thus, the tracking rate for Zq is a + 8. To estimate for the tracking rate for z~o, we need 
so called cone invariance property in [4]. It is a stronger version of the squeezing property 
(where the gap condition does not hold), which was originally introduced for the Navier- 
Stokes equation in [12] and improved in [8]. 
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Figure 4. (A) The long time behavior of (4.6) for two nearly differ- 
ent initial conditions, z\ = (0,0.02,0.01) andz 2 = (0,-0.02,0.01). (B) 
Leaves in the stable foliation through z\ and z 2 . (C) a\, b\, c\ £ ~$ Zl . 
(D) a 2 , b 2 , c 2 £ J( Z2 . 



Proposition 5.1. ( Cone Invariance Property) If zi ^ Z2 and denote u(t) = x{t,z\ ) — x(t 
andv(t) =y(t,z\)—y{t,Z2), then either 

(i) ||v(0||<Hf)||.V»GR 

(ii) \\v(t)\\ > \\u(t)\\,Vt£R 
( Hi) 3 t o £ R such that 




(b) (c) 



Figure 5. (A) For (4.6), z = (0.12,0.5,0.5) and fo = (0.12,0.5,0.0) 
and the tracking initial condition Zq = (0.2603,0.5225,-0.00364), 
found as the fixed point of £ 6 . (B) For (4.6). (C) For 8-mode KSE. 
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Lemma 5.2. For any zo € Z, either 

(5.2) \\z(t,zo)-z(t,z )\\ <^ a+5 >|M,V/>0, 
or there exists a to € (0,°°) such that 

(5.3) || Z ( M o)- Z (Mo)|| < | eW _ 5)f+(o(2W))||xo||) Vf e (foH _ 

Proof. Denote w(f) =x(/,Zo) — x(t,zo) and v(f) = y(f,Zo) — y(t>Zo)- It i s clear that 
0= ||v(0)| < ||m(0)|| = |jxo|j and thus by Proposition 5.1, either 

(i) ||v(t)||<H0ll, VreR 
or (ii) there exists to <G R such that 

'||v(f)|| < ||«(0ll , Vf <f 
|v(f)|| = |K0ll ,Vf = f 

|v(0ll > IKOH ,v/>/b. 

By the variational of constants, one has 



x(t,zo) =e ,A xo + [' e { '-^ A F{z{s,zo)) ds, andx(f,z~ ) = ( e^ A F{z{s,z )) ds. 
Jo Jo 

holds, then 

||z(f,Z0)-z(f,zb)|| = ||K(OII <e m \\x \\+S ['e^- s)A \\u(s)\\ds. 

Jo 



By the Gronwall inequality, one has that 

||«(0ll <e {a+5), \\xo\\,yt>0. 

This proves the first part of the lemma. 

Now if (ii) holds, for t < to, by the above estimate, we have 

\\z(t,zo)-z(t,z~o)\\<e ia+S)t \\*o\\- 
Using the Gronwall inequality and the fact that ||«(fo)| = ll v '( f o)|| we have for t > to, 

\\z(t,zo)-z(t,z~o)\\ = \\v(t)\\ < \\v(to)\\e(P- SH '-^ = \\u(to)\\e^^ 

< e («+ 5 )'||xo|| e ^- 5 »('-"' , = e « 3 -^' +r( '( 25+a -^||xo||. 

□ 

This lemma suggests that the tracking rate for fo is the same as the one for Zq over a 
short period of time, namely [0,fo] f° r some fo, but after that, the tracking rate becomes 
(J3 - 8)t +/o(25 + a - J3). It is easy to verify that 

(a + S)t < (P-S)t + t (2S + a-P), Vr > t Q . 

To show this result numerically, we consider both (4.6) with y = 25, dim(X) = 1 and 
dim(Y) = 2 with the initial condition zo chosen at random within a ball of radius 0.4 and the 
8-mode Galerkin approximation of the KSE with y = 32, dim(X) = 4, and dim(Y) = 4 and 
the initial condition zo chosen at random from a ball within radius 0.5. Figure 5(B) and (C) 
compare the rates of attractions for the projected and tracking initial data. Moreover, we 
compute their tracking rates via the simple linear regressions for these data sets in a short 
period of time (0.1 and .008 for (4.6) and 8-mode Galerkin approximation of the KSE) 
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(A) (4.6) (B) 8-mode KSE 



Figure 6. Tracking rate for tracking initial condition (x) is near the 
lower bound of the gap while the other (+) is near the upper bound of 
the gap. 

and they are -88.40 and 12.59 for (4.6) and -1640.45 and -400.04 for 8-mode Galerkin 
approximation of the KSE. The gap in eigenvalues is the interval (—99,21) for (4.6) and 
(—1700, —512) for 8-mode Galerkin approximation of the KSE. To ensure that this is not 
a special case, we repeat the above process with 10000 randomly chosen initial conditions 
and plot the tracking rates as shown in Figure 6. Observe that the data for tracking initial 
conditions are mainly close to the lower bound of the gap, while the data for projected 
initial conditions are mainly close to the upper bound of the gap. 
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